program test

call joke
!use mod_fftw

!  integer,parameter::M=15,N=15
!  real(8):: real2d(M,N)
!  complex(8):: complex2d(M/2+1,N)
!integer*8,save:: plan2dr2c,plan2dc2r
!
!
! ! call dfftw_init(M,N,real2d,complex2d)
!
!
!  call dfftw_plan_dft_c2r_2d(plan2dc2r,M,N,complex2d,real2d,FFTW_MEASURE)
!  do i=1,M/2+1
!  do j=1,N
!  complex2d(i,j)=cmplx(1D0,2D0+dble(j))
!  enddo
!  enddo
!
!  call dfftw_execute_dft_c2r(plan2dc2r,complex2d,real2d)
!  call dfftw_destroy_plan(plan2dc2r)
!  !call dfftw_c2r(real2d,complex2d)
!  !call dfftw_destroy()


endprogram

subroutine joke()
use imsl
use mod_upward
integer,parameter::MN=400
real(8)::mu,sigma,t,theta,pi,r
integer::i,j
real(8),dimension(:,:),allocatable::x,y,z,hz,s

allocate(x(MN,MN),y(MN,MN),z(MN,MN),hz(MN,MN),s(MN,MN))
pi=3.14159265358979
mu=4*pi*1D-7
sigma=1D0/100D0

t=1D-4
theta=sqrt(mu*sigma/4/t)

r=derf(0D0)

do i=1,MN
do j=1,MN
x(i,j)=(dble(i)-1-MN/2)*10D0
y(i,j)=(dble(j)-1-MN/2)*10D0
z(i,j)=100.D0*exp(-sqrt(x(i,j)**2+y(i,j)**2)/300.D0)
r=sqrt(x(i,j)**2+y(i,j)**2)+epsilon(1D0)
hz(i,j)=-(mu*sigma*r**3)**(-1)*(3*derf(theta*r)-2D0/sqrt(pi)*theta*r*(3+2*theta**2*r**2)*exp(-theta**2*r**2))
enddo
enddo

where(hz .gt. 1000D0) hz=0

where(hz .lt. -1000D0) hz=0

call init_upward(MN,MN,10.D0,10D0,z,r)

call timestamp()
do i=1,10
call equivsource(MN,MN,hz,s)
call timestamp()
enddo

endsubroutine

